function fsb_realign_multiple(file_names,rea_par)

% FSB - SPM:  realign multiple subjects/sessions in y-direction first and
% subsequently realign mean subject/session images over all directions
% Uses SPM2 for realignment
% Looks for the first underscore of the subject file name to determine the
% length of the name prefix string which it subsequently uses to decide
% which subject/session a volume belongs to
%
% EXAMPLE:
% fsb_realign_multiple(file_names,rea_par)
%
% INPUT
% A list of file names (file_names)
% A parameter struct rea_par with the fields
% rea_par.rp
% rea_par.srt
%
% rea_par.rp == 0:  Simple realigment
%                   Images are realigned in y-direction only, then a mean
%                   image is written for each session/run/trial without
%                   reslicing the other images first, so the mean images
%                   are not yet corrected for y-shifts.
%
%
% rea_par.rp > 0 :  Complex realignment
%                   Images are realigned in y-direction only, then a SPM
%                   Realignment is done and the spatial transformation
%                   combined. Mean images are therefore corrected for
%                   y-shifts.
%
% rea_par.rp == 1 : 2-step realignment and reslicing
%                   Mean images are realigned, the parameters for the
%                   y-direction realignment and the mean image realignment
%                   are combined on a volume-by volume basis and
%                   a mat file is saved for every volume separately
%                   Images are then resliced and written such that they can
%                   be loaded into fMRI Sandbox again
%
% rea_par.rp == 2 : Normalization and reslicing
%                   All images are normalized to the first volume in a
%                   series
%
% rea_par.rp == 3 : 2-step realignment and normalization
%                   Mean images are realigned, the parameters for the
%                   y-direction realignment are stored, and then the mean
%                   images are normalized to the first mean image as a
%                   template. The normalization parameters are saved and
%                   applied to all the volumes individually different.
%                   All images are normalized and written such that they
%                   can be loaded into fMRI Sandbox again.
%
% rea_par.rp == 4 : 2-step realignment and unwarping (not functional yet)
%                   Mean images are realigned, the parameters for the
%                   y-direction realignment and the mean image realignment
%                   are combined on a volume-by volume basis and
%                   a mat file is saved for every volume separately
%                   Images are then unwarped and written such that they
%                   can be loaded into fMRI Sandbox again.
%
% rea_par.srt:      Determines which column of the sandbox.intrial file is
%                   used for choosing the granularity of the realignment
%                   steps (e.g. trial-by-trial vs. session-by-session)
% Can use:
% A weighting image to determine which parts of the image should be used to
% determine the transformations.
%
% OUTPUT:
% Realigned images
% NOTES:
%
% External files:
% spm2_realign.m : Does the actual realignment in SPM2, just renamed
% fsb_load_files.m : Simple file loading function
% fsb_realing_yonly.m: Realignment in phase only
%
%
%$ Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin<2
    rea_par.rp = 3;
    rea_par.srt = 4;
end

set_spm2_path;
%set_spm5_path;
rehash;
global defaults
warning off; % Suppress warning about flipped image orientation

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Get files
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if ~exist('file_names','var');
    file_names = spm_get;
end

filnum = size(file_names,1);
file_names = cellstr(file_names);
file_names = file_names';
waitstr = ['Loading ' num2str(filnum) ' files'];
namecount = 1;
[pathstr, name{1}, ext, versn] = fileparts(file_names{1});
namepos = strfind(name,'_'); % Look for occurrences of underscores
namepos = namepos{1}; % Determine last letter of name prefix

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Determine file names
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

file_names_2{namecount}.name = name{1}(1:namepos); % Take the first namepos
%letters of the file name to determine to which dataset it belongs
file_names_2{namecount}.rname = ['r' name{1}(1:namepos) '.mat'];
file_names_2{namecount}.ons = 1;
file_names_2{namecount}.pathstr = pathstr;
thedir = file_names_2{1}.pathstr;
thedir = [thedir '\'];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Sort file prefixes
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

h = waitbar(0,waitstr);
for x = 1:filnum
    waitbar(x/filnum)
    [pathstr, name{x}, ext, versn] = fileparts(file_names{x});
    if x>1 & strcmp(name{x}(1:namepos),name{x-1}(1:namepos)) ==  0;
        namecount = namecount+1;
        file_names_2{namecount}.name = name{x}(1:namepos);
        file_names_2{namecount}.rname = ['r' name{x}(1:namepos) '.mat'];
        file_names_2{namecount}.ons = x;
        file_names_2{namecount}.pathstr = pathstr;
    end
end

close (h);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Run Realignment and Reslice in y-direction for each subject/session
% separately and get mean image names
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if rea_par.rp >0; % if this parameter is set to zero, mean images are directly produced

    meanfiles = {};
    for ii = 1:namecount;
        thename = file_names_2{ii}.name;
        all_files = dir([thedir thename '*.img']);
        all_files = vertcat(all_files(:).name);
        all_files = [all_files repmat(',1',size(all_files,1),1)];
        f = [repmat(thedir,size(all_files,1),1) all_files];
        if rea_par.rp ~= 2
            G = spm_figure('GetWin','Graphics');
            F = spm_figure('GetWin','Interactive');
            fg = spm_figure('FindWin','Interactive');
            spm_progress_bar('init');
            fsb_realign_yonly(f,0); %Change the second parameter to 1 if reslicing is necessary at this stage
            spm_progress_bar('clear');
            meanname = ['mean' file_names_2{ii}.name;];
        else
            meanname = ['meant' num2str(ii) '_'];
        end

        mean_files = dir([thedir meanname '*.img']);
        meanfiles{ii} = [thedir mean_files.name];
    end

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Get weighting images
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

meanfiles = char(meanfiles);
wimagename = [thedir 'wimage.img'];
wtest{1} = wimagename;

try
    fsb_load_files(wtest);
catch
    disp('no weighting image present, using standard procedures')
    wimagename = '';
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Run Realignment of mean images in all directions
% and normalize images to first image
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% wrap was [0 0 0] before

G = spm_figure('GetWin','Graphics');
F = spm_figure('GetWin','Interactive');
fg = spm_figure('FindWin','Interactive');
spm_progress_bar('init');

realign_flags = struct('quality',1,...
    'fwhm',2,...
    'sep',2,...
    'interp',2,...
    'wrap',[0 0 0],...
    'rtm',0,...
    'PW',wimagename,...
    'graphics',1); % set flags for realignment

% Determine realignment parameters and write mat files
P_rea = spm2_realign(meanfiles,realign_flags);

spm_progress_bar('clear');
warning off;
disp('done with mean image coregistration');

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Coregister all session images with realigned mean images of session
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

namecount = 1;

if rea_par.rp >0 && rea_par.rp ~=2;

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % If all images should be resliced before creating mean image
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    for x = 1:filnum

        if x == 1;
            mat_struct(x,:,:) = P_rea(namecount).mat;
        end

        if x>1

            if strcmp(name{x}(1:namepos),name{x-1}(1:namepos)) ==  0;
                namecount = namecount+1;
                %newname2 = [file_names_2{namecount}.pathstr '\' name{x+1} '.mat'];
                otto.mat(2,4) = P_rea(1).mat(2,4);
                newname2 = [file_names_2{namecount}.pathstr '\' name{x} '.mat'];
                %                 newname_mean = [meanfiles(namecount,1:end-4) '.mat'];
                %                 save (newname_mean,'mat');
            else
                newname2 = [file_names_2{namecount}.pathstr '\' name{x} '.mat'];
                otto = load(newname2,'mat');
            end

            % Save some parameters for later inspection
            newotto(x,:,:) = otto.mat;
            mat_struct(x,:,:) = P_rea(namecount).mat;

            % Calculate shift in y direction over scans
            startdev = otto.mat(2,4)-P_rea(1).mat(2,4);

            % Add shift to overall changes between sessions
            mat_struct(x,2,4) = P_rea(namecount).mat(2,4)+startdev;
            mat = squeeze(mat_struct(x,:,:));

            if rea_par.rp == 1 ||rea_par.rp == 5;
                % Change the mat files only if no normalization is done.
                % If Normalization is intended, just leave
                % the y-aligned mat files, and let normalization
                % do the rest
                save(newname2,'mat');
            end

        end
        newP(x,:,:) = P_rea(namecount).mat;
    end
else

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % If mean images should be created from non resliced data
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    matty = P_rea(2).mat;
    for x = 1:filnum

        if x>1 & strcmp(name{x}(1:namepos),name{x-1}(1:namepos)) ==  0;
            namecount = namecount+1;
        end

        newname = [file_names_2{namecount}.pathstr '\' name{x} '.mat'];
        mat = P_rea(namecount).mat;

        if newname>0
            save(newname,'mat');
        end
    end
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Prepare for reslicing or normalization
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
P = {};

for ii = 1:namecount;

    thename = file_names_2{ii}.name;
    all_files = dir([thedir 'r' thename '*.img']);

    if isempty(all_files) % if resliced images have not been written
        all_files = dir([thedir thename '*.img']);
    end

    all_files = vertcat(all_files(:).name);
    all_files = [all_files repmat(',1',size(all_files,1),1)];
    P{ii} = [repmat(thedir,size(all_files,1),1) all_files];
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Now either normalize all images to the mean of the first image
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if rea_par.rp == 3 || rea_par.rp ==2

    spm2_normalize(meanfiles,P,wimagename);

else
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Simply reslice all images
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    spm_progress_bar('clear');
    spm_progress_bar('init');
    reslice_flags = struct('interp',4,'mask',1,'mean',1,'which',2,'wrap',[1 1 0]'); % reslice images
    spm_reslice(P,reslice_flags);
    spm_progress_bar('clear');

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Try to unwarp images if so desired
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if rea_par.rp ==4

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Prepare for unwarping
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    if ~exist('matty','var')
        matty = squeeze(newP(2,:,:));
    end

    all_files1 = dir([thedir 'rt*.img']);
    all_files = vertcat(all_files1(:).name);
    P = [repmat(thedir,size(all_files,1),1) all_files];

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Define parameters for unwarping
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    par = struct('order',           [10 10],...
        'sfield',          [],...
        'M',               matty,...
        'sf_acq',          'Average',...
        'regorder',        1,...
        'lambda',          1e5,...
        'jm',              0,...
        'fot',             [2],...
        'sot',             [],...
        'fwhm',            4,...
        'rem',             1,...
        'exp_round',       'Average',...
        'noi',             5,...
        'hold',            [1 1 1 0 1 0]);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Try to run unwarping
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    try
        % Try to run unwarping
        ds = spm_uw_estimate(P,par);
        spm_uw_apply(ds);
    catch
        % Print message if not successful
        fprintf('Unwarping was not successful:\n%s',lasterr);
    end

end

close(fg);
close(G);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Plot the timecourses of the shifts
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
try
    if rea_par.rp>0
        figure(101);clf;
        subplot(3,1,1)
        plot(mat_struct(2:end,1,4));
        hold on;
        plot(newotto(2:end,1,4),'g');
        plot(newP(2:end,1,4),'r');
        subplot(3,1,2)
        plot(mat_struct(2:end,2,4));
        hold on;
        plot(newotto(2:end,2,4),'g');
        plot(newP(2:end,2,4),'r');
        subplot(3,1,3)
        plot(mat_struct(2:end,3,4));
        hold on;
        plot(newotto(2:end,3,4),'g');
        plot(newP(2:end,3,4),'r');
    end
end

disp('Realignment of mean images: Done');
end

